Admixture Proportions

All Samples

library(tidyverse)
library(fs)

ngsadmix_dir <- "data/ngsadmix_all/maf_0.05" # set NGSadmix outputs directory
N_K <- 10    # set number of K run
N_reps <- 4  # set number of reps run

# pull all log files
log_files <- list.files(ngsadmix_dir, pattern = ".log", full.names = T, recursive=T)

# read in all logs
all_logs <- lapply(1:length(log_files), FUN = function(i) readLines(log_files[i]))

# make list of the line that starts with "best like=" from all logs, just target 'b'
library(stringr)
bestlikes_str_list <- sapply(1:length(log_files), FUN= function(x) all_logs[[x]][which(str_sub(all_logs[[x]], 1, 1) == 'b')])

# make dataframe with 1:N_K and N_reps to add likelihood values
loglikes <- data.frame(K = rep(2:N_K, each=N_reps))

# add the log likelihood (first number in the string)
loglikes$loglike<-as.vector(as.numeric( sub("\\D*(\\d+).*", "\\1", bestlikes_str_list) ))

tapply(loglikes$loglike, loglikes$K, FUN= function(x) mean(abs(x))/sd(abs(x)))

           2            3            4            5            6            7            8            9           10 
1.314224e+03 3.724947e+08          Inf          Inf 5.753390e+02 6.269528e+03 1.961828e+04          Inf          Inf 

alewife_palette <- c("#636c62", "#869ca8", "#a2878a", "#b6bc9f", "#ecc2a3", "#eceeed", "#242b35", "#a39faa", "#ac9b7c")

metadata <- read_csv("data/herring_metadata.csv") %>%
  select(NMFS_DNA_ID, 
         GENUS, 
         SPECIES, 
         STATE_F, 
         WATERSHED, 
         WATER_NAME, 
         grouping_v3)

Rows: 94 Columns: 12── Column specification ───────────────────────────────────
Delimiter: ","
chr (11): sample, NMFS_DNA_ID, BOX_ID, BOX_POSITION, SA...
dbl  (1): BATCH_ID
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

gord <- c("BLUE", 
           "HYBR", 
           "GRTL", 
           "FINL", 
           "MIDL", 
           "CONL", 
           "NATLA", 
           "MIDA")
word <- c("Petitcodiac River", 
           "Hudson River", 
           "Lake Yonah", 
           "Lake Hartwell", 
           "Altamaha River", 
           "Roanoke River", 
           "Lake Superior", 
           "Lake Michigan", 
           "Lake Ontario", 
           "Canandaigua Lake", 
           "Cayuga Lake", 
           "Seneca Lake", 
           "Otisco Lake", 
           "East Grand Lake", 
           "Lake Champlain", 
           "Pattagansett Lake", 
           "Rogers Lake", 
           "Quonnipaug Lake", 
           "Miramichi River", 
           "Saco River", 
           "Black Creek", 
           "Choptank River")

labels <- read_tsv("data/admix-labels.tsv") 

Rows: 23 Columns: 4── Column specification ───────────────────────────────────
Delimiter: "\t"
chr (4): grouping_v3, group_label, WATER_NAME, pop_label
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

tmp <- metadata %>%
  mutate(xpos = 1:n())

group_pos <- tmp %>%
  group_by(grouping_v3) %>%
  summarise(midx = (min(xpos) - 0.5 + max(xpos) + 0.5)/ 2) %>%
  mutate(midy = 1)

group_labels <- group_pos %>%
  left_join(labels %>% group_by(grouping_v3) %>% slice(1))

Joining with `by = join_by(grouping_v3)`

K_2_rep_1 <- read.table("data/ngsadmix_all/maf_0.05/K_2_rep_1/output.qopt_with_sample_names", 
                        header = FALSE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

Error in `inner_join()`:
! Join columns in `x` must be present in the data.
✖ Problem with `sample`.
Backtrace:
 1. ... %>% arrange(gfact, wfact)
 5. dplyr:::inner_join.data.frame(., metadata, by = c(sample = "NMFS_DNA_ID"))

K_3_rep_1 <- read.table("data/ngsadmix_all/maf_0.05/K_3_rep_1/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_3_rep_2 <- read.table("data/ngsadmix_all/maf_0.05/K_3_rep_2/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_3_rep_3 <- read.table("data/ngsadmix_all/maf_0.05/K_3_rep_3/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_3_rep_4 <- read.table("data/ngsadmix_all/maf_0.05/K_3_rep_4/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

barplot(t(K_3_rep_1[2:4]), 
        col = alewife_palette, 
        names = K_2_rep_1$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 2, rep 1", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_3_rep_2[2:4]), 
        col = alewife_palette, 
        names = K_2_rep_2$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 2, rep 2", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_3_rep_3[2:4]), 
        col = alewife_palette, 
        names = K_2_rep_3$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 2, rep 3", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_3_rep_4[2:4]), 
        col = alewife_palette, 
        names = K_2_rep_4$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 2, rep 4", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

K_4_rep_1 <- read.table("data/ngsadmix_all/maf_0.05/K_4_rep_1/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_4_rep_2 <- read.table("data/ngsadmix_all/maf_0.05/K_4_rep_2/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_4_rep_3 <- read.table("data/ngsadmix_all/maf_0.05/K_4_rep_3/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_4_rep_4 <- read.table("data/ngsadmix_all/maf_0.05/K_4_rep_4/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

barplot(t(K_4_rep_1[2:5]), 
        col = alewife_palette, 
        names = K_4_rep_1$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 4, rep 1", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_4_rep_2[2:5]), 
        col = alewife_palette, 
        names = K_4_rep_2$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 4, rep 2", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_4_rep_3[2:5]), 
        col = alewife_palette, 
        names = K_4_rep_3$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 4, rep 3", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_4_rep_4[2:5]), 
        col = alewife_palette, 
        names = K_4_rep_4$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 4, rep 4", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

K_5_rep_1 <- read.table("data/ngsadmix_all/maf_0.05/K_5_rep_1/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_5_rep_2 <- read.table("data/ngsadmix_all/maf_0.05/K_5_rep_2/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_5_rep_3 <- read.table("data/ngsadmix_all/maf_0.05/K_5_rep_3/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_5_rep_4 <- read.table("data/ngsadmix_all/maf_0.05/K_5_rep_4/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

barplot(t(K_5_rep_1[2:6]), 
        col = alewife_palette, 
        names = K_5_rep_1$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 5, rep 1", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_5_rep_2[2:6]), 
        col = alewife_palette, 
        names = K_5_rep_2$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 5, rep 2", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_5_rep_3[2:6]), 
        col = alewife_palette, 
        names = K_5_rep_3$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 5, rep 3", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_5_rep_4[2:6]), 
        col = alewife_palette, 
        names = K_5_rep_4$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 5, rep 4", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

K_6_rep_1 <- read.table("data/ngsadmix_all/maf_0.05/K_6_rep_1/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_6_rep_2 <- read.table("data/ngsadmix_all/maf_0.05/K_6_rep_2/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_6_rep_3 <- read.table("data/ngsadmix_all/maf_0.05/K_6_rep_3/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_6_rep_4 <- read.table("data/ngsadmix_all/maf_0.05/K_6_rep_4/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

barplot(t(K_6_rep_1[2:7]), 
        col = alewife_palette, 
        names = K_6_rep_1$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 6, rep 1", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_6_rep_2[2:7]), 
        col = alewife_palette, 
        names = K_6_rep_2$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 6, rep 2", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_6_rep_3[2:7]), 
        col = alewife_palette, 
        names = K_6_rep_3$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 6, rep 3", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_6_rep_4[2:7]), 
        col = alewife_palette, 
        names = K_6_rep_4$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 6, rep 4", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

K_7_rep_1 <- read.table("data/ngsadmix_all/maf_0.05/K_7_rep_1/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_7_rep_2 <- read.table("data/ngsadmix_all/maf_0.05/K_7_rep_2/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_7_rep_3 <- read.table("data/ngsadmix_all/maf_0.05/K_7_rep_3/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_7_rep_4 <- read.table("data/ngsadmix_all/maf_0.05/K_7_rep_4/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

barplot(t(K_7_rep_1[2:8]), 
        col = alewife_palette, 
        names = K_7_rep_1$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 7, rep 1", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_7_rep_2[2:8]), 
        col = alewife_palette, 
        names = K_7_rep_2$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 7, rep 2", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_7_rep_3[2:8]), 
        col = alewife_palette, 
        names = K_7_rep_3$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 7, rep 3", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_7_rep_4[2:8]), 
        col = alewife_palette, 
        names = K_7_rep_4$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 7, rep 4", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

K_8_rep_1 <- read.table("data/ngsadmix_all/maf_0.05/K_8_rep_1/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_8_rep_2 <- read.table("data/ngsadmix_all/maf_0.05/K_8_rep_2/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_8_rep_3 <- read.table("data/ngsadmix_all/maf_0.05/K_8_rep_3/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_8_rep_4 <- read.table("data/ngsadmix_all/maf_0.05/K_8_rep_4/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

barplot(t(K_8_rep_1[2:9]), 
        col = alewife_palette, 
        names = K_8_rep_1$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 8, rep 1", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_8_rep_2[2:9]), 
        col = alewife_palette, 
        names = K_8_rep_2$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 8, rep 2", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_8_rep_3[2:9]), 
        col = alewife_palette, 
        names = K_8_rep_3$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 8, rep 3", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_8_rep_4[2:9]), 
        col = alewife_palette, 
        names = K_8_rep_4$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 8, rep 4", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

K_9_rep_1 <- read.table("data/ngsadmix_all/maf_0.05/K_9_rep_1/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_9_rep_2 <- read.table("data/ngsadmix_all/maf_0.05/K_9_rep_2/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_9_rep_3 <- read.table("data/ngsadmix_all/maf_0.05/K_9_rep_3/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_9_rep_4 <- read.table("data/ngsadmix_all/maf_0.05/K_9_rep_4/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

barplot(t(K_9_rep_1[2:10]), 
        col = alewife_palette, 
        names = K_9_rep_1$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 9, rep 1", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_9_rep_2[2:10]), 
        col = alewife_palette, 
        names = K_9_rep_2$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 9, rep 2", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_9_rep_3[2:10]), 
        col = alewife_palette, 
        names = K_9_rep_3$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 9, rep 3", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_9_rep_4[2:10]), 
        col = alewife_palette, 
        names = K_9_rep_4$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 9, rep 4", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

K_10_rep_1 <- read.table("data/ngsadmix_all/maf_0.05/K_10_rep_1/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_10_rep_2 <- read.table("data/ngsadmix_all/maf_0.05/K_10_rep_2/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_10_rep_3 <- read.table("data/ngsadmix_all/maf_0.05/K_10_rep_3/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

K_10_rep_4 <- read.table("data/ngsadmix_all/maf_0.05/K_10_rep_4/output.qopt_with_sample_names", 
                        header = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID")) %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)

barplot(t(K_10_rep_1[2:11]), 
        col = alewife_palette, 
        names = K_10_rep_1$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 10, rep 1", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_10_rep_2[2:11]), 
        col = alewife_palette, 
        names = K_10_rep_2$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 10, rep 2", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_10_rep_3[2:11]), 
        col = alewife_palette, 
        names = K_10_rep_3$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 10, rep 3", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

barplot(t(K_10_rep_4[2:11]), 
        col = alewife_palette, 
        names = K_10_rep_4$WATER_NAME, 
        las = 2, 
        cex.names = 0.5, 
        main = "K 10, rep 4", 
        space = 0, 
        border = "white", 
        cex.axis = 0.75, 
        cex.lab = 1, 
        xpd = NA)

ngsadmix_files <- dir_ls("data/ngsadmix_all", 
                         recurse = TRUE, 
                         glob = "*.qopt_with_sample_names")

ngsAdmix_tib <- lapply(ngsadmix_files, function(x) {
  read.table(x, 
             header = TRUE) %>%
    pivot_longer(cols = -sample, 
                 names_to = "Qval", 
                 values_to = "value") %>%
    mutate(path = x, 
           .before = sample)
}) %>%
  bind_rows() %>%
  filter(!is.na(value)) %>%
  mutate(Qval = str_replace(Qval, 
                            "X", 
                            "Q")) %>%
  extract(path, 
          into = c("K", 
                   "rep"), 
          regex = ".*K_([0-9]+)_rep_([0-9]+)/.*$", 
          convert = TRUE) %>%
  inner_join(., 
             metadata, 
             by = c("sample" = "NMFS_DNA_ID"), 
             relationship = "many-to-many") %>%
  mutate(gfact = factor(grouping_v3, 
                        levels = gord), 
         wfact = factor(WATER_NAME, 
                        levels = word)) %>%
  arrange(gfact, wfact)